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ABSTRACT 

High-energy gamma ray emission has been detected recently from supernovae remnants (SNRs) and their 
surroundings. The existence of molecular clouds near some of the SNRs suggests that the gamma rays originate 
predominantly from p-p interactions with cosmic rays accelerated at a closeby SNR shock wave. Here we 
investigate the acceleration of cosmic rays and the gamma ray production in the cloud self-consistently by 
taking into account the interactions of the streaming instability and the background turbulence both at the 
shock front and in the ensuing propagation to the clouds. We focus on the later evolution of SNRs, when the 
conventional treatment of the streaming instability is valid but the magnetic field is enhanced due to either 
Bell's current instability and/or due to the dynamo generation of magnetic field in the precursor region. We 
calculate the time dependence of the maximum energy of the accelerated particles. This result is then used 
to determine the diffusive flux of the runaway particles escaping the shock region, from which we obtain 
the gamma spectrum consistent with observations. Finally, we check the self-consistency of our results by 
comparing the required level of diffusion with the level of the streaming instability attainable in the presence 
of turbulence damping. The energy range of cosmic rays subject to the streaming instability is able to produce 
the observed energy spectrum of gamma rays. 

Subject headings: acceleration of particles-supernova remnants-instabilities-MHD-scattering-turbulence- 
cosmic rays. 



1. INTRODUCTION 

Cosmic rays below the knee are generally believed to 
originate from shock acceleration at the supernova remnants 
(SNRs). Although the process was proposed decades ago, it 
has not been easy to constrain the model from the observation- 
ally. Recent advances on gamma ray observations provide us 
a possibility to study the acceleration and ensuing propagation 
processes in detail. In particular, the molecular clouds in the 
vicinity of a SNR are good testbeds because of their high den- 
sity. Indeed, clouds illuminated by cosmic rays accelerated at 
a nearby SNR can be bright gamma ray sources (Aharonian 
& Atoyan 1996; Fatuzzo et al. 2006; Gabici et al. 2009). 

Gamma rays are produced through the decay of neutral pi- 
ons generated in p-p collisions (hadronic origin) or through 
inverse Compton interactions by energetic electrons with am- 
bient photon target fields such as the cosmic microwave 
background radiation (leptonic origin). It is still under de- 
bate whether the observed gamma ray emissions are of the 
hadronic or leptonic origin although the hadronic scenario is 
favored by many authors because of the evidence of strongly 
amplified magnetic fields near SNRs (see the review by Blasi 
2010 and references therein). 

We consider the gamma ray emission from massive molec- 
ular clouds near SNRs, motivated, in particular, by the re- 
cently detected gamma ray emission from the SNR W28. As 
an old SNR, the emission from W28 is perceived of purely 
hadronic origin since the life time of high energy electrons 
is short. The observed gamma emission from W28 indicates 
that the flux of cosmic rays is much enhanced there compared 
to the typical Galactic values. This places a constraint on the 
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cosmic ray diffusion. The enhancement of scattering corre- 
sponding to a reduction (1 — 10%) of the CR diffusion co- 
efficient compared to Galactic mean has been suggested by 
a few earlier models to match the observations (Fujita et al. 
2009, Ohira et al. 2011, Li & Chen 2010). No physical 
justification on how and why the scattering is boosted has 
been provided, however. In addition, earlier studies made a 
few serious assumptions including the Bohm diffusion and a 
phenomenological power law evolution of the maximum en- 
ergy accelerated at a certain epoch. We shall reinvestigate 
this problem by incorporating a proper physical description 
of the relevant processes, i.e. the shock acceleration in the 
presence of the streaming instability and nonlinear damping 
processes by background turbulence. The slower diffusion 
implies enhanced wave perturbations, which may arise from 
streaming instability (see, e.g., Longair 2002). Earlier work 
has shown that streaming instability is limited by background 
turbulence (Yan & Lazarian 2002; Farmer & Goldreich 2004; 
Beresnyak & Lazarian 2008). This could be the determinative 
factor for the maximum energy attainable at the shock front 
as suggested by Ptuskin & Zirakashvili (2005). However, the 
isotropic Kolmogorov scaling for the turbulence, adopted by 
Ptuskin & Zirakashvili (2005), is not applicable to MHD tur- 
bulence. Moreover, if the enhanced scattering in the vicinity 
of SNRs as indicated by the observations are due to the in- 
creased flux of cosmic rays there, we need to examine whether 
the flux of accelerated particles is sufficient to induce high 
enough growth rates of the streaming instability to overcome 
the nonlinear damping by background turbulence. 

In this paper, we apply our present day understanding of 
the interaction between the streaming instability and the back- 
ground turbulence to the modeling of the gamma ray emission 
from molecular clouds near SNRs. We shall treat the problem 
in a self-consistent way by comparing the streaming level that 
is allowed by the preexisting turbulence and the required dif- 
fusion for the CRs. In §2 we determine the maximum energy 
of accelerated particles at shocks, in §3 we calculate the spec- 
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trum of cosmic rays in the vicinity of SNRs, in §4 we obtain 
the spectrum of gamma ray produced by the p-p interactions 
of the escaping cosmic rays in a nearby cloud and compare 
it with observations, §5, we check the self-consistency of our 
result by comparing the required the diffusion near SNRs and 
the streaming level allowed in the presence of turbulence. The 
discussion and summary are provided in §6 and §7. In the Ap- 
pendix we provide a list of notations used in the paper. 

2. MAXIMUM ENERGY OF CRS ACCELERATED AT THE SHOCKS 

Diffusive shock acceleration of energetic CR particles relies 
on the crucial process of amplification of MHD turbulence so 
that particles can be trapped at the shock front long enough to 
be accelerated to the high energy observed. One of the most 
popular scenarios that has been adopted in the literature is the 
streaming instability generated by the accelerated particles. 
However, in the highly nonlinear regime the fluctuations of 
magnetic field arising from the streaming instability get large 
and the classical treatment of the streaming instability is not 
applicable. We circumvent the problem by proposing that the 
field amplification we consider does not arise from the stream- 
ing instability, but is achieved earlier through other processes, 
e.g. the interaction of the shock precursor with density per- 
turbations preexisting in the interstellar medium (Beresnyak, 
Jones, & Lazarian 2009). Due to the resonant nature of the 
streaming instability, the perturbations SB arising from it are 
more efficient in scattering CRs compared to the large scale 
fluctuations produced by non-resonant mechanisms, e.g. the 
one in Beresnyak et al. (2009). Therefore in this paper, we 
limit our discussions to the regime of bB < Bq, where Bq is 
the magnetic field that has already been amplified in the pre- 
cursor region 4 . 

When particles reach the maximum energy at a certain time, 
they escape and the growth of the streaming instability stops. 
Therefore we can obtain the maximum energy by considering 
the stationary state of the evolution. The steady state energy 
density of the turbulence W(k) at the shock is determined by 



(t/±v A )W(fe) =2{T cr -T d )W{kl 



(1) 



where U is the shock speed, and the term on the l.h.s. rep- 
resents the advection of turbulence by the shock flow. v A = 
Bo/y4nnm and n are the Alfven speed and the ionized gas 
number density of the precursor region, respectively. The plus 
sign represents the forward propagating Alfven waves and the 
minus sign refers to the backward propagating Alfven waves. 
The terms on the r.h.s. describes the wave amplification by 
the streaming instability and damping with r fr ,r^ as the cor- 
responding growth and damping rates of the wave. The distri- 
bution of accelerated particles at strong shocks is f(p) p~ A . 
If taking into account the modification of the shock structure 
by the accelerated particles, the CR spectrum becomes harder. 
Assume the distribution of CRs at the shock is fo(p) °^ p~ 4+a . 
The nonlinear growth was studied by Ptuskin & Zirakashvili 
(2005). 

It was demonstrated by Schlickeiser & Shalchi (2008) that 
waves can grow or damp at the shock precursor depending 
on the the spatial boundary conditions. If assuming an equal 
amount of forward and backward waves at the shock front, the 
forward wave will be growing and the backward wave will 
be efficiently damped in the upstream region. We therefore 

4 The effective Bo is therefore renormalized and can be much larger than 
the typical field in ISM (see, e.g., Diamond & Malkov 2007). 



neglect the backward moving wave mode and consider here 
only the growing forward moving mode, which is the one that 
effectively contributes to the particle scattering. The general- 
ized growth rate of streaming instability then is 



Un 2 q 2 v A Vl+A 2 
c^k 



x / dpp 



D 



dx 



(2) 



where q is the charge of the particle, c is the light speed, 
p res = ZeBoVl+A 2 /c/k res is the momentum of particles that 
resonate with the waves. A = hB/Bo is wave amplitude nor- 
malized by the mean magnetic field strength Bq. 



D = y/l +A 2 vr g /3/A 2 (> k res ) (3) 

is the diffusion coefficient of CRs, v and r g are the velocity 
and Larmor radius of the CRs. The distribution function of 
CRs is 
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where \ measures the ratio of CR pressure at the shock and 
the upstream momentum flux entering the shock front, m is 
the proton rest mass, and p max is the maximum momentum 
accelerated at the shock front. H(p) is the Heaviside step 
function. 

In the planar shock approximation, the distribution of ac- 
celerated particles is 



/i(p,*)=/o(p)exp( -(U 



V A, 



D(p,xi) 



(5) 



at the upstream of the shock (x > 0) and / = /o at the down- 
stream. Insert Eqs.(l3l l4l5l l into Eq.(f2]l, one gets the following 
growth rate of the upstream forward moving wave at x=0, 



r cr (fe) = 



C c £U 2 (U + v A )k 



i-o 



(l+A2)(l-«)/2cVA<|>(j> B 



(6) 



where C cr = 9/2/(4 - a)/ (2 - a), r = mc 2 /q/B . The lin- 
ear damping is negligible since the medium should be highly 
ionized. In fully ionized gas, there is nonlinear Landau 
damping, which, however, is suppressed due to the reduc- 
tion of particles' mean free path in the turbulent medium (see 
Yan & Lararian 2011; Brunetti & Lazarian 2011); we there- 
fore neglect this process here. Background turbulence itself 
can cause nonlinear damping to the waves (Yan & Lazarian 
2002). Unlike hydrodynamical turbulence, MHD turbulence 
is anisotropic with eddies elongated along the magnetic field. 
The anisotropy increases with the decrease of the scale (Gol- 
dreich & Sridhar 1995). At the scales of the Larmor radii Tl 
of the cosmic rays, which are also the wavelengths 1 jk of the 
waves induced by the streaming instability, the scale disparity 
becomes very large with k± ^ fe| , where we use fej| to distin- 
guish the parallel component of the turbulence wave packet 
wavenumber fej from the parallel wavenumber of the growing 

wave feii, In MHD turbulence, wave packet cascades after it 
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travels a distance of ~ L 1 < /3 £ ± 2 ' /3 ^> On the other 

hand, the instability grows fastest for the most parallel wave 
(k ~ &||) allowed in a turbulence medium with their wave num- 
bers satisfying k±/k^ ~ 8B/B ~ (fcxL) -1 / 3 . Because of the 
scale disparity, > &x ^ &|, the nonlinear damping rate in 
MHD turbulence is less than the wave frequency k^v A , and it 
is given by (Farmer & Goldreich 2004; Yan & Lazarian 2004; 
Beresnyak & Lazarian 2008) 

r d ~y/k/Lv A , (7) 

where L is the injection scale of background turbulence, and 
the k is set by the resonance condition k ~ k\\ ~ 1 / r^. Inserting 

Eqs.lgE) into Eq.Q and adopting UdW /dx « f/ 2 W/Z) in the 
case of efficient wave amplification, one gets 



3U 2 A 2 
2v(l +A 2 ) 



C„^C/ 2 (f/ + v A ) 



cvA^(p„ m )(^o)«(l+A 2 )( 1 -")/ 2 ' 



(8) 

There are various models for the diffusive shock acceler- 
ation. We consider here the escape-limited acceleration. In 
this model, particles are confined in the region near the shock 
where turbulence is generated. Once they propagate far up- 
stream at a distance I from the shock front, where the self- 
generated turbulence by CRs fades away, the particles escape 
and the acceleration ceases. The characteristic length that 
particles penetrate into the upstream is D(p)/U. The maxi- 
mum momentum is reached when D{p)/U ~ 1/4 5 . Assuming 
I oc R sht then the maximum momentum of particles acceler- 
ated during the Sedov phase is determined by the condition 



D(p max )=KUR sh , 



(9) 



wher e K < 1 is a numerical factor, see table 2. From equations 
J3I9K we get 

p max 3KA 2 UR sh 



mc VT +A 2 vr Q 
Insert Eq.([T0jl into Eq.®, we get for A < l 6 



(10) 



A = 



where §(pmax) is approximated by (p m ax/mc) a /a. 
In the limit of low shock velocity, 



v A < U < c 



V.A 



1 1/4 
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we get 



3\2 



2 4 ' 

r c l v A 



(12) 



(13) 



5 The factor 1/4 arises from the following reason. As pointed out by Os- 
trowski & Schlickeiser (1996), the spectrum is steepened for small 1, i.e., 
lU/D(p) < 4 

6 We neglect a factor of \/l +A 2 here. 



10' 



Energy of CRs release at time t after SN explosion 



-Our results 
-Simple power law 
Earlier result 




FIG. 1 . — The energy of CRs that are released at the shock at time t in the 
Sedov phase. Our result shows that the often assumed power law solution (see 
Gabici et al. 2009; Ohiraet al. 2010) is only realized in asymptotic regime 
as described in Eqs. fll2H"3"V It is also larger than the earlier result (dotted 
line) in Ptuskin & Zirakashvili (2005) where the damping of the waves by 
background turbulence is overestimated. 

At the Sedov phase (f > t sed = 250(£ , 5 i/(n t/ 9 5 )) 1 ' /3 yr), where 
£51 = £'sAr/10 51 erg and Ug = C/// 10 9 cm/s are the total energy 
of explosion and the initial shock velocity, the evolution of 
shock radius and speed are governed by 



R sh =4.3(e 51 /« ) 1/5 f^. 5 P c , 
[7 = 1.7 x 10 3 (e 5 i/«o) 1/5 ^ /5 km/s. 



(14) 




From Eqs. dl3ll4b . we see that Pmax °^ t > and 8 — 18/5. In 
FigQ] we plot pmax/ {mc) vs. the time t since supernova ex- 
plosion. The solid line represents the results from Eqs.dTTl. 
As we see, at earlier epoch when advection and streaming in- 
stability are both important, the evolution of p max does not 
follow a power law. For comparison, we also put a power law 
evolution in the same figure as depicted by Eq.([T3l (dashed 
line). Our result is also larger than that obtained by Ptuskin & 
Zirakashvili (2005) since the wave dissipation rate is overes- 
timated in their treatment. 
2 

3. THE COSMIC RAY DISTRIBUTION NEAR SNRS 

> At a given time, only particles with p > p max can escape the 
shock front. Adopting the simplified approximation, 



< RshU, p<p„ 
> R sh U, p > p„ 



(15) 



at a given energy E = cp, there is a one to one correspon- 
dence between the CR momentum and the CR escape time 
t esp (or radius R esp = (1 + K)R S h) in spite of the fact that the 
acceleration is continuous during the shock expansion since 
Pmaxit) reaches p. If the maximum momentum has a power 



law dependence on t, one can easily gets t e 



-1/8 



Since we consider the later stage of SNR acceleration, the 
shock radius cannot be neglected and the point source as- 
sumption is not applicable in general. One needs to take into 
account the spatial distribution of the sources. Nevertheless, 
in the case that the diffusion distance R c i = 2^D(t age — t esp ) 
is larger than the radius R esp , the distribution function of en- 
ergetic proton at a distance r can be described as if the CRs 
were from a point source (see, e.g., Ohira et al. 201 1) 
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F(E) 



f(E) 



exp 



(16) 



The diffusion coeffcient D = %Dism, where Djsm = 
10 28 v /17T0GeVcm 2 /s (Gabici et al. 2009; Berezinskii et al. 
1990). f(E) is the spectrum of runaway particles that is inte- 
grated over shock expansion. 

We adopt the routine from Ohira et al. (2010) here. If the 
accelerated particles at a given time is 



f(p,t)dpdt = K(t) 



( p \ - 2+a 


P 


1 — exp 


Pmaxif) _ 


Vmc/ 



d\ntdp, 



(17) 

then the general spectrum of protons dispersed from the ac- 
celerator is 



f(p)< 



p- 1+ "K( Pm Up)) 
( {p)[d Pnmx /dt} t=pmUp) 



(18) 



where p m l x {p) is the inverse function of p mil x(t)- K(t) is a 
normalization factor and can be estimated from 



4ft rPmaAt) 

P C R~^ dpvp 3 f (p) = tymll 2 . 



(19) 



In the case of flatter spectrum at the shock front, i.e., a > 0, 
K(t) oc U 2 R 3 sh p~". Evaluation of dp max /dt from Eq.dTTb 
shows that it can be well represented by a power law in spite of 
the fact that p max does not follow an exact power law. Insert 
K(t) and dp max /dt ~ Pmax/t into Eq[l8] we get a universal 
power spectrum f(E) oc t\EsnE~ s with s = 2 regardless of the 
flatter original spectrum at the shock front 7 . r| is the fraction 
of SN energy converted into CRs. We obtain a similar cos- 
mic ray momentum spectrum as earlier authors (e.g., Ptuskin 
& Zirakashvili 2005; Ohira et al. 2010) without assuming the 
power law evolution of p max (,t). Pmax{t) in our paper was ob- 
tained by considering the balance of the wave growth from the 
streaming instability, the advection by the shock flow as well 
as wave damping by the background turbulence. The latter is 
due to the nonlinear interaction between the wave and the tur- 
bulence, and has not been properly accounted for. Ptuskin & 
Zirakashvili (2005) considered a similar process, but assumed 
isotropy of the MHD turbulence which is not valid. 

Although time evolution of p max {t) does not follow a power 
law at high energies 8 , the final CR spectrum F(E) at a 
later age of Sedov phase and a distance r > Rj is well re- 
produced by the power law form if the diffusion coefficient 
D(p) has a power law dependence on particle momentum 
p. This is because high energy CRs are produced at earlier 
epoch t esp <C t age so that R c [ can be well approximated by 
R c i ~ 2^jDt uge . Adopting the set of parameters for W28 as 
listed in Tables 1,2, we get the flux of CRs E 2 F(E). In FigH 
we plot the flux of CRs at different epochs after the SN ex- 
plosion. All the values of the parameters except that for tur- 
bulence are compatible with those from earlier work ( Fujita 
et al. 2009; Gabici et al. 2010; Li & Chen 2010). We adopt 
a turbulence injection scale of 3Qpc, which is consistent with 
the scenario that interstellar turbulence is injected from super- 
novae explosions. The CR energy at each peak corresponds 

7 The influence of CR cooling on the spectrum index is negligible (Ohira 
et al. 2010). 

8 It is a good approximation for E max < lOOOGeV (or t > 4kyr) in the case 
we consider (see Fig[T}. 
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FIG. 2. — The spectrum of CRs at a distance r = 12pc after 1800 (solid 
line), 6000 (dotted line), 12000 (dashdot line), 50000 years (cross line). The 
Galactic mean is plotted as a reference (dashed line). 

to the maximum momentum p max (t) reached at each epoch t. 
The galactic background is also plotted for comparison. Ap- 
parently, the CR intensities are dramatically enhanced in the 
vicinity of SNRs as compared to the Galactic mean intensi- 
ties. 

4. GAMMA RAY FLUX IN THE VICINITY OF SNRS 

In the case F(E) follows a power law distribution, the pion 
gamma ray emissivity is given by (see Aharonian & Atoyan 
1996) 



q(E y ) 



-OppF(Ey)a\A, 



(20) 



where Ey « iip/10 is the corresponding yray energy, G pp ps 
30 x [0.95 + 0.06 \n(E p/GeV )]mb is the cross section for pp 
collisions at E p , f K is the fraction of energy that is transferred 
from parent protons to secondary pions, r|^ ~ 1 .4 — 1 .5 is a 
parameter to account for the contribution from both cosmic 
rays and the interstellar gas (Dermer 1986). The total flux 
then is 



dNy _ M c q(E) 



dE y 



4nd 2 r, 



(21) 



where M c , d are the mass and the distance of the cloud. Com- 
bining Eqs.( |16|2QI2U . we obtain the flux of gamma ray emis- 
sion as shown in FigfJ) where our result is plotted against both 
Fermi and H.E.S.S. data. The GeV and TeV data are adopted 
from Abdo et al. (2010) and Aharonian et al. (2008), respec- 
tively. Our result produces a power law spectral index ~ 2.75, 
showing that the steepening of the spectrum of particle can be 
naturally explained by the propagation effect. Indeed, similar 
fits have been also obtained by other models, e.g., Ohira et 
al (2011), Li & Chen (2010). We did not need, however, ei- 
ther to assume Bohm diffusion and phenomenological power 
law evolution of the momentum of the escaping partciels as 
in Ohira et al. (201 1), or a steeper spectrum for the escaping 
particles as in Li & Chen (2010). A decrease of the spatial 
diffusion coefficient by a factor of % = 0.05 comparable to 
earlier works (Gabici 2010, Li & Chen 2010) is inferred here, 
which can originate from streaming instability. Estimates will 
be provided below to demonstrate that the instability can grow 
for the energy range we consider. 

5. ENHANCED SCATTERING AND STREAMING INSTABILITY 
NEAR SNRS 
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TABLE 1 
Physical parameters of W28 
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FIG. 3. — The spectrum of Gamma ray emission from W28. The Fermi 
data are shown as dotted points (Abdo et al. 2010), and the H.E.S.S. data are 
plotted as 'x' points (Aharonian et al. 2008) with error bars. Solid line is our 
result. 

Our results show that the local scattering of CRs has to be 
enhanced by an order of magnitude % = 0.05 in order to pro- 
duce the amount of y ray emission observed. A natural way 
to increase the scattering rate is through the streaming insta- 
bility. We provide here a self-consistency check by exam- 
ining whether the streaming instability operates in the pres- 
ence of nonlinear damping by the background turbulence. We 
envisage a local low density, high temperature cavity with 
Bear = 2/jG surrounding the supernova subjected to strong UV 
radiation and stellar wind (see Fugita et al. 2009), so that there 
is no ion-neutral damping in this case. Only nonlinear damp- 
ing exists (see eq[7j. The growth rate in the linear regime is 

9 



N(>E) 



VA 



(22) 



where v s is the streaming speed of CRs. The growth rate 
should overcome the damping rate (eqj7j for the instability 
to operate. The condition Y gr > F c / leads to 



v s >v A \ 1 



hva 



(23) 



The spatial diffusion coefficient adopted here, D s» v s L = 
%Djsm, satisfies this requirement. The growth and damping 
rates are compared in Figj4] We see that the streaming in- 
stability works in the energy range needed to produce the 
observed y ray emission, proving that our results are self- 
consistent. 

Note that the case we consider here is different from the 
general interstellar medium discussed in Yan & Lazarian 
(2004) and Farmer & Goldreich (2004), namely, the local cos- 
mic ray flux near SNRs is much enhanced (see Fig. 2). Conse- 
quently, the growth rate of the streaming instability becomes 
high enough to overcome the damping rate by the preexisting 
turbulence in the considered energy range. 

9 We neglect the nonlinear Landau damping, which is suppressed in turbu- 
lence due to decrease of mean free path. 
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FIG. 4. — The growth and nonlinear damping rates of streaming instability. 
With the locally enhanced flux, the growth rate of streaming instability be- 
comes much larger than the mean Galactic value so that it can overcome the 
nonlinear damping by turbulence for a wide energy range. This is consistent 
with our earlier treatment in which streaming instability plays an essential 
role in the cosmic ray diffusion near SNRs. 



TABLE 2 

Model parameters adopted in the paper 



1 



K 



a 



0.1-0.3 -0.05 -0.3 0.04-0.1 0.2-0.4 0.5 



6. DISCUSSION 

Gamma rays have been detected from molecular clouds 
near SNRs. They are believed to the result of interactions with 
the energetic particles accelerated at SNR shocks, and thus 
provide the best probe to both the shock acceleration accel- 
eration and propagation processes of cosmic rays near SNRs. 
To determine the spectrum of CRs that generates the gamma 
ray emissions, it is necessary to establish the time dependence 
of the energy of escaping particles. Earlier it has been done 
assuming that only CRs with p > p max {t) can escape and that 
the time dependence is a power law (see Gabici et al 2010, 
Ohira et al. 2011). In this paper, we derive our results on 
the basis of the well motivated and tested physical model of 
the streaming instability in the presence of background turbu- 
lence. With this model we obtain p max (t) which is determined 
by the interaction of streaming instability at the upstream and 
the background turbulence. We limit the discussions to the 
later evolution of SNR acceleration which is stage of the ob- 
served gamma ray sources. In this case, the classical formula 
for the growth of streaming is valid. 

The spectral fit at high gamma ray energies with H.E.S.S. 
data is less optimal. High energy cosmic rays are generated 
earlier when strong magnetic amplification is needed. As we 
pointed out earlier, the classical treatment of streaming is in- 
valid and moreover, the field amplification may well be due to 
other processes, eg., the current instability (Bell 2004) and/or 
turbulent density perturbations interacting with the shock pre- 
cursor (Beresnyak et al. 2009). 

Another crucial ingredient that determines the CR flux at 
the clouds is the spatial diffusion coefficient of CRs escap- 
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ing from SNRs. It has been speculated recently that this spa- 
tial diffusion is suppressed due to CR induced instabilities to 
match the observed level of gamma ray fluxes (see Gabici et 
al 2010). No specific study has been provided so far. Here we 
place a constraint to the level of streaming by applying our 
present understanding of turbulence damping of waves (Yan 
& Lazarian 2004, Farmer & Goldreich 2004, Beresnyak & 
Lazarian 2008). By balancing the growth rate of the stream- 
ing instability with the turbulence damping rate, we obtain 
not only the energy range at which the streaming instability 
operates, but also the lower limit of the CR spatial diffusion 
coefficient. These results are found to be consistent with the 
required enhanced scattering indicated from the interpretation 
of the gamma ray data. 

For the spectral index of the the CR power law spectrum, 
we took into account the hardening of the accelerate parti- 
cles at the shock front. This does not contradict the fact that 
the spectrum of runaway particles is steeper as it is a time- 
averaged result. For the same reason, we do not consider that 
the decrease of the compression ratio would contribute to the 
steepening of the power law spectrum of the accelertaed CRs. 
In fact, as long as s < 2, the spectral index of runaway parti- 
cles becomes uniformly 2. The main steepening at a distance 
r from the SNR is due to the propagation effect. In the case 
when r>R c [> R esp , the point source is a good approximation 

for SNR, and the spectrum is steepened by p 3a / 2 if D °< p a . 

The diffusion coefficient D near the shock is assumed to 
be of the form D = %10 2S ^/E/T0GeVcm 2 /s following ear- 
lier treatments (see Gabici et al 2010, Ohira et al. 2011). 
Undoubtedly, solving the diffusion coefficient surrounding 
SNRs taking into account both streaming instability and back- 
ground turbulence is an important step toward a complete self- 
consistent picture. It will be one of our future endeavors. 
Another simplification we made concerns the magnetic field 
strength at the shock. For the later stage of Sedov phase evolu- 
tion, we assumed an already increased magnetic field strength 
<~ 100/jG, and did not consider the specifics of earlier ampli- 
fication process, which is widely believed to be present. This 
process is important for earlier shock evolution and acceler- 
ation of ultra high energy CRs. It, however, is a subject of 
intensive debates and is definitely beyond the scope of our 
current paper. 

The set of physical parameters and model parameters we 
used are listed, respectively, in tables 1&2. The physical con- 
dition we adopted for W28 is close to that used in literatures 
( Fujita et al. 2009; Gabici et al. 2010; Li & Chen 2010). 
As for the model parameters, we have tight constraints only 
for the diffusion coefficient (a, %), similar to those obtained in 
earlier literatures (Gabici et al. 2010; Li & Chen 2010) and 



that for the energy of the SN explosion r|. The results are not 
so sensitive to the other model parameters, as shown in table 
2. 

7. SUMMARY 

We investigated the acceleration of particles at the later 
stage of SNR evolution and the escape of these particles. We 
calculated the flux of the escaping CRs at a molecular cloud 
in the vicinity of SNRs. We quantitatively took into account 
the competition of turbulence generation by the streaming in- 
stability of CRs and background turbulence damping, both, 
near the shock waves and in the ensuing propagation of CRs 
from SNRs. The resulting gamma ray spectrum from the p-p 
interactions is compared with both Fermi and H.E.S.S. data. 
Our main results are: 

• The streaming instability plays a crucial role for cos- 
mic ray acceleration, particularly at the later stage of 
SNRs. The competition of turbulence generation by 
the streaming instability and turbulence damping deter- 
mines the maximum energy of the accelerated particles. 

• The spectrum of runaway particles follows a universal 
power law of f(E) °= E~ 2 if the original spectrum index 
s at the shock front is s < 2 for the later Sedov phase. 
The main steepening at a distance from the SNR is due 
to the propagation effect. 

• The flux of CRs is increased by several orders of mag- 
nitude compared to the mean Galactic value, creating 
enough turbulence by the streaming instability in the 
vicinity of the shock which overcomes the damping by 
background turbulence for the energy range considered. 

• The enhanced scattering by the instability is enough to 
reproduce the observed gamma ray flux from massive 
molecular clouds near the SNR. 
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A 


Normalized wave amplitude SB/Bq 


a 


hardening of the CR spectium at the shock front 


Bo 


mean magnetic field at the shock in the later Sedov phase 


Bcav 


inercloud magnetic field strength 


SB 


wave amplitude 


c 


light speed 


d 


distance of the molecular cloud from observer 


D 


diffusion coefficient of CRs 


E 


CR energy 




supernova explosion 


f 


distribution function of CRs 


fx 


fraction of energy transferred from parent protons to pions 


k 


wave number 


K(t) 


Normalization factor of CR distribution function 


L 


the injection scale of background turbulence 


m 


proton rest mass 


M c 


cloud mass 


n 


intercloud number density 


Ny 


y ray flux 


P 


CR's momentum 




the maximum momentum accelerated at the shock front 


PCR 


CR pressure 


q 


charge of the particle 


r 


distance from SNR centre 


R c 


the distance of the molecular cloud from the SNR centre 


r g 


Larmor radius of CRs 


Rd 


diffusion distance of CRs 


R s h 


shock radius 




the escaping distance/time of CRs 


s 


ID spectrum index of CR distribution 


t 


time since supernova explosion 


f 


the age of SNR 


t7ed 


the time at which SNR enters the Sedov phase 


U 


shock speed 


Ui 


initial shock velocity 


V 


particle speed 


Vj 


streaming speed of CRs 


W 


wave energy 


a 


power index of D with respect to particle momentum p 


X 


reduction factor of D with respect to D/sm 


8 


power index of p max with respect to t 


n 


fraction of SN energy converted into CRs 




a numerical factor in Eq|20l 


r„ 


the growth rate of streaming instability 


r rf 


wave damping rate 


K 


ratio of diffusion length to shock radius 


Qo 


the Larmor frequency of nonrelativistic protons 




cross section for pp collision 


k 


the ratio of CR pressure to fluid ram pressure 



APPENDIX 
LIST OF NOTATIONS 

The notation we used is listed in Table lAl 



